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Abstract. We review the development of update schemes for quantum lattice models simulated us- 
ing world line quantum Monte Carlo algorithms. Starting from the Suzuki-Trotter mapping we dis- 
cuss limitations of local update algorithms and highlight the main developments beyond Metropolis- 
style local updates: the development of cluster algorithms, their generalization to continuous time, 
the worm and directed-loop algorithms and finally a generalization of the flat histogram method of 
Wang and Landau to quantum systems. 



QUANTUM MONTE CARLO WORLD LINE ALGORITHMS 

Suzuki's realization in 1976 [1] that the partition function of a c/-dimensional quantum 
spin-1/2 system can be mapped onto that of a (d+ 1) -dimensional classical Ising model 
with special interactions enabled the straightforward simulation of arbitrary quantum 
lattice models, overcoming the restrictions of Handscomb's method [2]. Quantum spins 
get mapped onto classical world lines and the Metropolis algorithm [3] can be employed 
to perform local updates of the configurations. 

Just like classical algorithms the local update quantum Monte Carlo algorithm suffers 
from the problem of critical slowing down at second order phase transitions and the 
problem of tunneling out of metastable states at first order phase transitions. Here we 
review the development of non-local update algorithms, stepping beyond local update 
Metropolis schemes: 

• 1993: the loop algorithm [4], a generalization of the classical cluster algorithms to 
quantum systems allows efficient simulations at second order phase transitions. 

• 1996: continuous time versions of the loop algorithm [5] and the local update al- 
gorithms [6] remove the need for an extrapolation in the discrete time step of the 
original algorithms (an approximation-free power- series scheme had been intro- 
duced for the S=l/2 Heisenberg model already in [2], and a related, more general 
method with local updates was presented in [7]). 

• from 1998: the worm algorithm [8], the loop-operator [9, 10] and the directed 
loop algorithms [11] remove the requirement of spin-inversion or particle-hole 
symmetry. 

• 2003: flat histogram methods for quantum systems [12] allow efficient tunneling 
between metastable states at first order phase transitions. 
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WORLD LINES AND LOCAL UPDATE ALGORITHMS 



The Suzuki-Trotter decomposition 

In classical simulations the Boltzmann weight of a configuration c at an inverse 
temperature /3 = l/ksT is easily calculated from its energy E c as exp(— fiE c ). Hence 
the thermal average of a quantity A 

(A)ciassicai = £A C exp(-/3£ c )/ £exp(-/3£ c ) (1) 

c c 

can be directly estimated in a Monte Carlo simulation. The key problem for a quantum 
Monte Carlo simulation is that the simple exponentials of energies get replaced by 
exponentials of the Hamilton operator if: 

(A) =Tr[Aexp(-j8H)]/Tr[exp(-j8H)] (2) 

The seminal idea of Suzuki [1], using a generalization of Trotter's formula [13], was 
to split H into two or more terms H = £f Hi so that the exponentials of each of the terms 
exp(— flHi) is easy to calculate. Although the //, do not commute, the error in estimating 
the exponential 

exp(-eH) ^T7_exp( -£//■) + €?(e 2 ) (3) 

is small for small prefactors e and better formulas of arbitrarily high order can be derived 
[14]. Applying this approximation to the partition function we get Suzuki's famous 
mapping, here shown for the simplest case of two terms Hi and H2 

Z = Tr[exp(-j6//)] = Tr [exp(-AT(//i +// 2 )] M 

= Tr[exp(-AT//i)exp(-AT// 2 )] M + €?(At 2 ) (4) 
= £ (h\Ui\i 2 )(i2\U2\i3)---(i2M-i\Ui\i m )(i2M\U 2 \h) + ^(At 2 ), 

'1i-i'2M 

where the time step is At = /3/M, the each are complete orthonormal sets of basis 
states, and the transfer matrices are £/, = exp(— At//,). The evaluation of the matrix 
elements {i\U\\i') is straightforward since the H t are chosen to be easily diagonalized. 



The World Line Representation 

As an example we consider a one-dimensional chain with nearest neighbor interac- 
tions. The Hamiltonian H is split into odd and even bonds H\ and H2, as shown in Fig. 
la). Since the bond terms in each of these sums commute, the calculation of the expo- 
nential is easy. Equation (4) can be interpreted as an evolution in imaginary time (inverse 
temperature) of the state by the "time evolution" operators U\ and U2. Within each 
time interval At the operators U\ and and U2 are each applied once. This leads to the fa- 
mous "checkerboard decomposition", a graphical representation of the sum on a square 
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FIGURE 1. The "checkerboard decomposition": a) the Hamiltonian is split into odd and even bond 
terms, b) A graphical representation of Suzuki's mapping of a one-dimensional quantum system to a 
two-dimensional classical one, where an example world line configuration is shown. 



lattice, where the applications of the operators U\ are marked by shaded squares (see Fig. 
lb). The configuration along each time slice corresponds to one of the states |4) in the 
sum (4). 

This establishes the mapping of a one-dimensional quantum to a two-dimensional 
classical model where the four classical states at the corners of each plaquette interact 
with a four-site Ising-like interaction. For Hamiltonians with particle number (or magne- 
tization) conservation we can take the mapping one step further. Since the conservation 
law applies locally on each shaded plaquette, particles on neighboring time slices can 
be connected and we get a representation of the configuration {\ik}} in terms of world 
lines. The sum over all configurations with non-zero weights (ijt|t/|ijk+i) corre- 

sponds to the sum over all possible world line configurations. In Fig. lb) we show such 
a world line configuration for a model with one type of particle (e.g. a spin- 1/2, hardcore 
boson or spinless fermion model). For models with more types of particles there will be 
more kinds of world lines representing different particles (e.g. spin-up and spin-down 
fermions). 



Local Updates 

The world line representation can be used as a starting point of a quantum Monte 
Carlo algorithm [15]. Since particle number conservation prohibits the breaking of world 
lines, the local updates need to move world lines instead of just changing local states as 
in a classical model. 

As an example we consider a one-dimensional tight binding model with Hamiltonian 

H=-tJ^(c}c i+1 +c] +1 c?) , (5) 
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FIGURE 2. Examples of the two types of local moves used to update the world line configuration in 
a tight-binding model with two states per site and Hamiltonian Eq. (5): a) plaquette weights (/jt|£/|ijt+i) 
of the six possible local world line configurations in a tight binding model; b) the two types of updates in 
discrete time and c) in continuous time. 



where c] creates a particle (spinless fermion or hardcore boson) at site i. Fig. 2a shows 
the plaquette weights (ijfc|E/|ijt+i) for each of the six world line configurations on a 
shaded plaquette in this model. 

The local updates are quite simple and move a world line across a white plaquette 
[15, 16], as shown in Fig. 2b). Slightly more complicated local moves are needed for 
higher-dimensional models [17], t-J models [18, 19] and Kondo lattice models [19]. 

Since these local updates cannot change global properties, such as the number of 
world lines or their spatial winding, they need to be complemented with global updates 
if the grandcanonical ensemble should be simulated [17]. The problem of exponentially 
low acceptance rate of such moves was remedied only much later by the non-local update 
algorithms discussed below. 



The Continuous Time Limit 

The systematic error arising from the finite time step At was originally controlled by 
an extrapolation to the continuous time limit At — > from simulations with different 
values of the time step At. It required a fresh look at quantum Monte Carlo algorithms 
by a Russian group [6] in 1996 to realize that, for a discrete quantum lattice model, this 
limit can already be taken during the construction of the algorithm and simulations can 
be performed directly at At — > 0, corresponding to an infinite Trotter number M = oo. 

In this limit the Suzuki-Trotter formula Eq. (4) becomes equivalent to a time- 
dependent perturbation theory in imaginary time [6, 8]: 



Trexp(-JBH) =Tr 



exp(-/3# )^exp / dxV(x) 
Jo 



(6) 



Tr 



exp(-j6tf ) - J q P dxV(x)dx + ^J* dx x jf P dx 2 V (ti)V(t 2 ) + . 



where the symbol & denotes time-ordering of the exponential. The Hamilto- 
nian H = Ho + V is split into a diagonal term Ho and an offdiagonal pertur- 
bation V. The time-dependent perturbation in the interaction representation is 
V(t) = exp(xHo)V exp(— xHq). In the case of the tight-binding model the hopping 
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term t is part of the perturbation V, while additional diagonal potential or interaction 
terms would be a part of Hq. 

To implement a continuous time algorithm the first change in the algorithm is to 
keep only a list of times at which the configuration changes instead of storing the 
configuration at each of the 2M time slices in the limit M — > <*>. Since the probability 
for a jump of a world line [see Fig. 2a)] and hence a change of the local configuration is 
sinh(Ar?) « Ar« 1 /M the number of such changes remains finite in the limit M — > <*>. 
The representation is thus well defined, and, equivalently, in Eq. (6) only a finite number 
of terms contributes in a finite system. 

The second change concerns the updates, since the probability for the insertion of a 
pair of jumps in the world line [the upper move in Fig. 2b)] vanishes as 

insert jump = sinh 2 (AT0/cosh 2 (AT0 - At 2 - 1/M 2 - (7) 

in the continuous time limit. To counter this vanishing probability, one proposes to insert 
a pair of jumps not at a specific location but anywhere inside a finite time interval [6]. 
The integrated probability then remains finite in the limit At — > 0. Similarly instead of 
shifting a jump by At [the lower move in Figs. 2b,c)] we move it by a finite time interval 
in the continuous time algorithm. 



Stochastic Series Expansion 

An alternative Monte Carlo algorithm, which also does not suffer from time dis- 
cretization, is the stochastic series expansion (SSE) algorithm [7], a generalization of 
Handscomb's algorithm [2] for the Heisenberg model. It starts from a Taylor expansion 
of the partition function in orders of /3 : 

°° Qn 

Z = Trexp(-/3tf) = £ ^ J Tr(-H) n 

<x> on 

= E^T E E (h\-H bl \i 2 )(i 2 \-H b2 \i 3 )...(i n \-H bn \h) (8) 

n=0 H - {i u ...i n }{b u ...b n } 

where in the second line we decomposed the Hamiltonian H into a sum of single-bond 
terms H = and again inserted complete sets of basis states. We end up with 

a similar representation as Eq. (4) and a related world-line picture with very similar 
update schemes. For more details of the SSE method we refer to the contribution of 
A.W. Sandvik in this proceedings volume. 

The SSE representation can be formally related to the world line representation by 
observing that Eq. (8) is obtained from Eq. (6) by setting Hq — 0, V = H and integrating 
over all times (compare also Fig. 3) T,- [20]. This mapping also shows the advantages 
and disadvantages of the two representations. The SSE representation corresponds to a 
perturbation expansion in all terms of the Hamiltonian, whereas world line algorithms 
treat the diagonal terms in Hq exactly and perturb only in the offdiagonal terms V of 
the Hamiltonian. World line algorithms hence need only fewer terms in the expansion, 
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FIGURE 3. A comparison of a) world lines in discrete time, b) in continuous time and c) a similar 
configuration in the SSE representation. In the SSE representation the continuous time index is replaced 
by an integer order index of the operators, at the cost of additional diagonal terms (the dashed lines). 



but pay for it by having to deal with imaginary times T ; . The SSE representation is thus 
preferred except for models with large diagonal terms (e.g. bosonic Hubbard models) or 
for models with time-dependent actions (e.g. dissipative quantum systems [21]). 



THE LOOP ALGORITHM 

While the local update world line and SSE algorithms enable the simulation of quantum 
systems they suffer from critical slowing down at second order phase transitions. Even 
worse, changing the spatial and temporal winding numbers has an exponentially small 
acceptance rate. While the restriction to zero spatial winding can be viewed as a bound- 
ary effect, changing the temporal winding number and thus the magnetization or particle 
number is essential for simulations in the grand canonical ensemble. 

The solution to these problems came with the loop algorithm [4] and its continuous 
time version [5]. These algorithms, generalizations of the classical cluster algorithms 
[22] to quantum systems, not only solve the problem of critical slowing down, but also 
updates the winding numbers efficiently for those systems to which it can be applied. 

Since there is an extensive recent review of the loop algorithm [23], we will only 
mention the main idea behind the loop algorithm here. In the classical Swendsen- 
Wang cluster algorithm each bond in the lattice is considered, and with a probability 
depending on the local configuration two neighboring spins are either "connected" or 
left "disconnected", as shown in Fig. 4a). "Connected" spins form a cluster and must 
be flipped together. Since the average extent of these cluster is just the correlation 



a) O 



1b) | f <p^> 



FIGURE 4. a) in the cluster algorithms for classical spins two sites can either be connected (thick line) 
or disconnected (thin line), b) in the loop algorithm for quantum spins two or fours spins on a shaded 
plaquette must be connected. 



Non-local Updates for Quantum Monte Carlo Simulations 



February 2, 2008 



6 



a) 



b) 



c) 



FIGURE 5. A loop cluster update: a) world line configuration before the update, where the world line 
of a particle (or up-spin in a magnetic model) is drawn as a thick line and that of a hole (down-spin) as a 
thin line; b) world line configuration and a loop cluster (grey line); c) the world line configurations after 
all spins along the loop have been flipped. 



length of the system, updates are performed on physically relevant length scales and 
autocorrelation times are substantially reduced. 

Upon applying the same idea to world lines in QMC we have to take into account 
that (in systems with particle number or magnetization conservation) the world lines 
may not be broken. This implies that a single spin on a plaquette cannot be flipped 
by itself, but at least two, or all four spins must be flipped in order to create valid 
updates of the world line configurations. Instead of the two possibilities "connected" or 
"disconnected", four connections are possible on a plaquette, as shown in Fig. 4b): either 
horizontal neighbors, vertical neighbors, diagonal neighbors or all four spins might be 
flipped together. The specific choices and probabilities depend, like in the classical 
algorithm, on details of the model and the world line configuration. Since each spin 
is connected to two (or four) other spins, the cluster has a loop-like shape (or a set of 
connected loops), which is the origin of the name "loop algorithm" and is illustrated in 
Fig. 5. 

While the loop algorithm was originally developed only for six- vertex and spin- 1/2 
models [4] it has been generalized to higher spin models [24], anisotropic spin models 
[25], Hubbard [26] and t-J models [27]. 



Applications of the loop algorithm 

Out of the large number of applications of the loop algorithm we want to mention only 
a few which highlight the advances made possible by the development of this algorithm 
and refer to Ref. [23] for a more complete overview. 

• The first application of the discrete and continuous time loop algorithms [28, 5] 
were high accuracy simulations of the ground state parameteres of the square lattice 
Heisenberg antiferromagnet, establishing beyond any doubt the existence of Neel 
order even for spin S = 1/2. 

• The exponential divergence of the correlation length in the same system could be 
studied on much larger systems with up to one million spins [29, 30, 31] and with 
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much higher accuracy than in previous simulations [17], investigating not only the 
leading exponential behavior but also higher order corrections. 

• For quantum phase transitions in two-dimensional quantum Heisenberg antiferro- 
magnets, simulations using local updates had been restricted to small systems with 
up to 200 spins at not too low temperatures and had given contradicting results 
regarding the universality class of the phase transitions [32, 33]. The loop algo- 
rithm enabled simulations on up to one hundred times larger systems at ten times 
lower temperatures, allowing the accurate determination of the critical behavior at 
quantum phase transitions [34, 35]. 

• Similarly, in the two-dimensional quantum XY model the loop algorithm allowed 
accurate simulations of the Kosterlitz-Thouless phase transition [36], again improv- 
ing on results obtained using local updates [37]. 

• In SU(4) square lattice antiferromagnets, the loop algorithm could clarify that a 
spin liquid state thought to be present based on data obtained using local update 
algorithms on small lattices [38] is actually Neel ordered [39]. 

• A generalization, which allows to study infinite systems in the absence of long 
range order, was invented [40]. 

• The meron cluster algorithm, an algorithm based on the loop algorithm, solves the 
negative sign problem in some special systems [41]. 

WORM AND DIRECTED LOOP ALGORITHMS 

Problems of the loop algorithm in a magnetic field 

As successful as the loop algorithm is, it is restricted - as the classical cluster algo- 
rithms - to models with spin inversion symmetry (or particle- hole symmetry). Terms in 
the Hamiltonian which break this spin-inversion symmetry - such as a magnetic field 
in a spin model or a chemical potential in a particle model - are not taken into account 
during loop construction. Instead they enter through the acceptance rate of the loop flip, 
which can be exponentially small at low temperatures. 

As an example consider two 5=1/2 quantum spins in a magnetic field: 

// = /SiS 2 -M^+>^) (9) 

In a field h = 7 the singlet state l/\/2(| fj) — | ||)) with energy —3/47 is degenerate 
with the triplet state | ft) with energy 1/47 — h = —3/47, but he loop algorithm is 
exponentially inefficient at low temperatures. As illustrated in Fig. 6a), we start from the 
triplet state | ft) an d propose a loop shown in Fig. 6b). The loop construction rules, 
which do not take into account the magnetic field, propose to flip one of the spins 
and go to the intermediate configuration | fj,) with energy —1/47 shown in Fig. 6c). 
This move costs potential energy 7/2 and thus has an exponentially small acceptance 
rate exp(— /37/2). Once we accept this move, immediately many small loops are built, 
exchanging the spins on the two sites, and gaining exchange energy 7/2 by going to the 
spin singlet state. A typical world line configuration for the singlet is shown in Fig. 6d). 
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FIGURE 6. A loop update for two antiferromagnetically coupled spins in a magnetic field with / = h. a) 
Starting from the triplet configuration 1 b) a loop is constructed, proposing to go to c), the intermediate 
configuration | ||), which has an exponentially small acceptance rate, and finally into configurations like 
d) which represent the sing let state 1/V2(| U) - I IT}). As in the previous figure a thick line denotes an 
up-spin and a thin line a down-spin. 



The reverse move has the same exponentially small probability, since the probability to 
reach a world line configuration without any exchange term [Fig. 6c)] from a spin singlet 
configuration [Fig. 6d)] is exponentially small. 

This example clearly illustrates the reason for the exponential slowdown: in a first 
step we lose all potential energy, before gaining it back in exchange energy. A faster 
algorithm could thus be built if, instead of doing the trade in one big step, we could 
trade potential with exchange energy in small pieces, which is exactly what the worm 
algorithm does. 



The Worm Algorithm 

The worm algorithm [8] works in an extended configuration space, where in addition 
to closed world line configurations one open world line fragment (the "worm") is 
allowed. Formally this is done by adding a source term to the Hamiltonian which for 
a spin model is 

tfworm = #-T7£(S++Sr)- (10) 

i 

This source term allows world lines to be broken with a matrix element proportional to 
r\. The worm algorithm now proceeds as follows: a worm (i.e. a world line fragment) is 
created by inserting a pair Sj) of operators at nearby times, as shown in Fig. 7a,b). 
The ends of this worm are then moved randomly in space and time [Fig. 7c)], using local 
Metropolis or heat bath updates until the two ends of the worm meet again as in Fig. 
7d). Then an update which removes the worm is proposed, and if accepted we are back 
in a configuration with closed world lines only, as shown in Fig. 7e). This algorithm 
is straightforward, consisting just of local updates of the worm ends in the extended 
configuration space but it can perform nonlocal changes. A worm end can wind around 
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FIGURE 7. A worm update for two antiferromagnetically coupled spins in a magnetic field with J = h. 
a) starting from the triplet configuration | ||) a worm is constructed in b) by inserting a pair of S + and S~ 
operators, c) these "worm end" operators are then moved by local updates until d) they meet again, when 
a move to remove them is proposed, which leads to the closed world line configuration e). As in the two 
previous figures a thick line denotes an up-spin and a thin line a down-spin. 



the lattice in the temporal or spatial direction and that way change the magnetization and 
winding number. 

In contrast to the loop algorithm in a magnetic field, where the trade between potential 
and exchange energy is done by first losing all of the potential energy, before gaining 
back the exchange energy, the worm algorithm performs this trade in small pieces, 
never suffering from an exponentially small acceptance probability. While not being 
as efficient as the loop algorithm in zero magnetic field (the worm movement follows 
a random walk while the loop algorithm can be interpreted as a self-avoiding random 
walk), the big advantage of the worm algorithm is that it remains efficient in the presence 
of a magnetic field. 

A similar algorithm was already proposed more than a decade earlier [42]. Instead of 
a random walk using fulfilling detailed balance at every move of the worm head in this 
earlier algorithm just performed a random walk. The a posteriori acceptance rates are 
then often very small and the algorithm is not efficient, just as the small acceptance rates 
for loop updates in magnetic fields make the loop algorithm inefficient. This highlights 
the importance of having the cluster-building rules of a non-local update algorithm 
closely tied to the physics of the problem. 



The Directed Loop Algorithm 

Algorithms with a similar basic idea are the operator-loop update [9, 10] in the SSE 
formulation and the directed-loop algorithms [11] which can be formulated in both 
an SSE and a world-line representation. Like the worm algorithm, these algorithms 
create two world line discontinuities, and move them around by local updates. The 
main difference to the worm algorithm is that here these movements do not follow an 
unbiased random walk but have a preferred direction, always trying to move away from 
the last change. The directed loop algorithms might thus be more efficient than the worm 
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algorithm but no direct comparison has been performed so far. For more details see the 
contribution of A.W. Sandvik in this volume. 



Applications 



Just as the loop algorithm enabled a break-through in the simulation of quantum mag- 
nets in zero magnetic field, the worm and directed loop algorithms allowed simulations 
of bosonic systems with better efficiency and accuracy. A few examples include: 

• Simulations of quantum phase transitions in soft-core bosonic systems, both for 
uniform models [8] and in magnetic traps [43]. 

• By being able to simulate substantially larger latttices than by local updates [44] 
the existence of supersolids in hard-core boson models was clarified [45] and 
the ground-state [45, 46] and finite-temperature phase diagrams [47] of two- 
dimensional hard-core boson models have been determined. 

• Magnetization curves of quantum magnets have been calculated [48]. 



FLAT HISTOGRAMS AND FIRST ORDER PHASE TRANSITIONS 



The main problem during the simulation of a first order phase transition is the exponen- 
tially slow tunneling time between the two coexisting phases. For classical simulations 
the multi-canonical algorithm [49] and recently the Wang-Landau algorithm [50] eases 
this tunneling by reweighting configurations such as to achieve a "flat histogram" in 
energy space. In a canonical simulation the probability of visiting an energy level E is 
p(E)p(E) oc p(£')exp(— fiE) where the density of states p(E) is the number of states 
with energy E. While the multi-canonical algorithm [49] changes the canonical distri- 
bution p(E) by reweighting it in an energy-dependent way, the algorithm by Wang and 
Landau discards the notion of temperature and directly uses the density of states to set 
p(E) \/p(E), which gives a constant probability in energy space p(E)p(E) = const.. 
The unknown quantity p(E) is determined self-consistently in an iterative way and then 
allows to directly calculate the free energy 



and other thermodynamic quantities at any temperature. The main change to a simulation 
program using a canonical distribution is to replace the canonical probability p(E) = 
exp(— j5E) by the inverse density of states p(E) = l/p(E). 

This algorithm cannot be straightforwardly used for quantum systems, since the 
density of states p(E) is not directly accessible for those. Instead we recently proposed 
[12] to start from the SSE formulation of the partition function Eq. (8): 



F = -fc s rin£p(£)exp(-/3£) 



(ID 



E 



n 



F = — 



& B rinTrexp(-j6#) = -k B T\n Y ^-Tr(-H) 



n 



n=0 
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= -k B T\n^— £ E (h\-H bl \i 2 )(i2\-H b2 \i 3 )---(i„\-H bn \i l ) 

n=0 H - {h,...i„} {bi,...b„} 

OO 

= -k B T\nY,P n g(n). (12) 

n=0 

The coefficient g(n) of the n-th order term in an expansion in the inverse temperature /3 
now plays the role of the density of states p(E) in the classical algorithm. Similar to the 
classical algorithm, by using 1 / g(n) as the probability of a configuration instead of the 
usual SSE weight, a flat histogram in the order n of the series is achieved. Alternatively 
instead of such a high-temperature expansion a finite-temperature perturbation series 
can be formulated [12]. 

This algorithm was shown to be effective at first order phase transitions in quantum 
systems and promises to be effective also for the simulation of quantum spin glasses. 



WHICH ALGORITHM IS THE BEST? 

Since there is no "best algorithm" suitable for all problems we conclude with a guide on 
how to pick the best algorithm for a particular problem. 

• For models with particle-hole or spin-inversion symmetry a loop algorithm is 
optimal [4, 5, 9]. Usually an SSE representation [9] will be preferred unless the 
action is time-dependent (such as long-range in time interactions in a dissipative 
quantum system) or there are large diagonal terms, in which case a world line 
representation is better. 

• For models without particle hole symmetry a worm or directed-loop algorithm is 
the best choice: 

- if the Hamiltonian is diagonally dominant use a worm [8] or directed loop [11] 
algorithm in a world line representation. 

- otherwise ause directed-loop algorithm in an SSE representation. [9, 10, 11]. 

• At first order phase transition a generalization of Wang-Landau sampling to quan- 
tum systems should be used [12]. 

The source code for some of these algorithms is available on the Internet. Sandvik 
has published a FORTRAN version of an SSE algorithm for quantum magnets [51]. 
The ALPS (Algorithms and Libaries for Physics Simulations) project is an open-source 
effort to provide libraries and application frameworks for classical and quantum lattice 
models as well as C++ implementations of the loop, worm and directed-loop algorithms 
[52]. 
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